Mathematica for Quantum Information Theory

The basics

This notebook was inspired by Jony Hudson’s Mathematica course (which is now rather hard to find online!)

Section 1: Pauli matrices and single qubit gates

Like we did in week 1, let’s define the Pauli matrices as

Id ={{1, 0}, {0, 1}};
X = {{0, 1}, {1, 0}};
Y = {{0, -I}, {I, 0}};
Z = {{1, 0}, {0, -1}};

paulis = {X, Y, Z};

We naughtily defined these as upper-case letters, which you should generally avoid doing so as not to collide with inbuilt functions like N. We could have done this much more concisely using the PauliMatrix function:

paulis = {X,Y,Z} = PauliMatrix[Range[3]];

To select the kth entry of a list X in Mathematica we can use X[[k]]. Note the indices start at 1 (not 0)! To print matrices nicely in Mathematica we can then use the MatrixForm function. e.g.

paulis[[2]] // MatrixForm

where // is just the postfix syntax for applying a function:

a // f

f @ a

Note that “directly” multiplying matrices...

X Y

X * Y

will multiply them element-wise, i.e. where (Subscript[(XY), ij]) = (Subscript[X, ij])(Subscript[Y, ij]). Matrix multiplication as we ordinarily define it is performed using Dot:

X.Y //MatrixForm

Let’s define a function to compute the commutation relation between two matrices.

commutator[A_, B_]:= A.B - B.A

commutator[X,Y] // MatrixForm

The operator := is SetDelayed, which you should typically always use (instead of just =) when defining functions. The latter evaluations the function definition immediately which often isn’t what we want!

Define a function that computes the anti-commutator and produce a Table (or take a more sophisticated approach if you like) to verify that the standard anti-commutation rules hold

Solution

Define a function that implements evolution under a Pauli operator by some angle. That is, create a function that evaluates (e^(-i θ P)) where θ is a scalar and P is either X, Y or Z

Hints:
1) Use MatrixExp.
2) Mathematica has a long list of helpful short cuts for symbols
e.g. see https://reference.wolfram.com/language/guide/MathematicalTypesetting.html and https://reference.wolfram.com/language/guide/SpecialCharacters.html

Solution

Generalize this to implement evolution about an arbitrary axis r on the Bloch Sphere, i.e. (e^(-i θ n .σ)) where n is any 3D vector and σ = {X, Y, Z}

Solution

Simplify the expression for (e^(-i θ n .σ)) when n is a real unit-vector.

Solution

Assessed problem

Define a function to create an arbitrary single qubit state and hence convince yourself that evolution under (e^(i α (OverHat[n]).σ)) implements a rotation by an angle -2α about unit vector (OverHat[n]) on the Bloch sphere.

Hint: Use the representation of the Bloch vector in spherical coordinates

Solution

Section 2: Multi-qubit systems and gates

Define projectors onto the 0 and 1 states of a qubit. Use these, and the KroneckerProduct function, to define the CNOT gate (controlling upon the left qubit)

Solution

The Hadamard gate is defined as follows

H = (1/(Sqrt[2])){{1,1},{1,-1}};
H // MatrixForm

It is used to prepare the plus and minus states, e.g.

H.{1,0}
H.{0,1}

Show that the Hadamard gates and CNOT gates can be combined to prepare a Bell state. (How would you prepare the other Bell states?)

Solution

Below, we define calcPartialTrace which accepts a matrix and a list of qubits (indexed from zero).

getPartialKet[bitInds_List, bitVals_List, numAllBits_] := Module[
	{ops = ConstantArray[IdentityMatrix[2], numAllBits]},
	Do[
		ops[[ bitInds[[i]]+1 ]] = If[bitVals[[i]]===0, {{1},{0}}, {{0},{1}}],
		{i, Length@bitInds}];
	KroneckerProduct @@ Reverse[ops]]

calcPartialTrace[ρ_?MatrixQ, qubits_List] /; Length@qubits < Log2@Length@ρ :=
	Sum[
		With[
			{ket = getPartialKet[qubits, vals, Log2@Length@ρ]},
			ConjugateTranspose[ket] . ρ . ket],
		{vals, Tuples[{0,1}, Length@qubits]}]

It can be used like so:

ρm = Table[ρ[i,j],{i,4},{j,4}];
ρm // MatrixForm

calcPartialTrace[ρm, {0}] // MatrixForm

χ = Table[RandomComplex[],(2^3),(2^3)];
ρm = χ . (χ^\[ConjugateTranspose])/Tr[χ . (χ^\[ConjugateTranspose])] // Chop;
ρm // MatrixForm

calcPartialTrace[ρm, {0,2}] // Chop // MatrixForm

Show that tracing out half of any Bell state generates the maximally mixed state

Solution

Assessed problem

Suppose you generate a generic entangled state with the following circuit:


Write a function to calculate the reduced state of each qubit given an arbitrary angle θ and unit vector n. Does it matter which qubit you look at?

Solution

Use Manipulate and Plot to study how the purity of the reduced state depends on the rotation angle ϑ  and the rotation direction r. Comment on what you observe. 
(Hint: think about the best way to represent r)

Solution